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Auxiliary field quantum Monte Carlo methods for Hubbard models are generally based on a 
Hubbard-Stratonovitch transformation where the field couples to the ^-component of the spin. This 
transformation breaks SU(2) spin invariance. The symmetry is restored only after summation over 
the auxiliary fields. Here, we analyze an alternative decomposition, which conserves SU(2) spin 
invariance, but requires the use of complex numbers. We show that this algorithm gets rid of the 
very large fluctuations observed in imaginary time displaced correlation functions of quantities which 
do not commute with the ^-component of the total spin. The algorithm prooves to be efficient for 
the study of spin dynamics. 

Auxiliary field quantum Monte Carlo Algorithms for Hubbard models are usually based on the discrete Hubbard- 
Stratonovtich decomposition |l| : 

cxp (-ArL^^T-fJ [ n ll~^)J =C 1^ exp \ a l^ s i{ n iA - n ll) I ' W 

Here, r\q a — ct cj CT where ct (c* CT ) creates (annihilates) an electron on site i with z-component of spin a, cosh(a) = 

exp(Arf7/2) On an TV-site lattice, the constant C = exp (AtUN/A) /2 and At corresponds to an imaginary time 
step. As apparent from the above equation, for a fixed set of Hubbard- Stratonovitch (HS) fields, Si . . . sn, SU (2)-spin 
symmetry is broken, (i.e. the expression is not invariant under the transformation cj a — > [exp (i<f>ea /2)] a s , cj g , . Here, 
e is a unit vector and a is a vector consisting of the Pauli-spin matrices.) Clearly SU(2) spin symmetry is restored 
after summation over the HS fields. 
Alternatively, one may consider [hj 




^-atu ( n ?, T - 5) ( n u - 1) j = c E ± ex p | w E s ? { n l 



rxp I -Arf'Y^ ( n fi -_)(„, -_)]= f y I ,s ; (/ l . T + n U -lJJ . (2) 

where cos(a) = exp (— AtII/2) and C = exp (AtUN/A) /2 N . With this choice of the HS transformation SU(2) spin 
invariance is retained for any given HS configuration. 

The aim of this note, is to address the question: will we obtain a more efficient and/or reliable quantum Monte Carlo 
algorithm if we enhance the number of symmetries conserved by the HS transformation. We consider the extended 
Hubbard model: 

i i i 

with the hopping kinetic energy 

K i=J2{ C h C Ul<r + 4 + l^)- (4) 

<7,S 

Here W > 0, S = ±& x , ±a y where a x , a y are the lattice constants. We impose twisted boundary conditions: 

C ?+La x , CT = 6X P (2™*/$o) % a , C I+LSy a = C la , (5) 

with $0 = hc/e the flux quanta and L the linear length of the square lattice. The boundary conditions given by Eq. 
(|s|) account for a magnetic flux threading a torus on which the lattice is wrapped. At half-filling, and constant value 
oiU/t the W-tevTH drives the ground state from a antiferromagnetic Mott insulator to a d x 2_ y 2 superconductor ||. 
At U/t = 4, this quantum transition occurs at W c /t ~ 0.3. At finite values of W < W c , numerical simulations are 
consistent with the occurrence of a d x 2_ y 2 superconductor upon doping of the Mott insulating state ||. 
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To decompose the W-term, we use the approximate relation 



\ J2 7(0exp(VA7w 7 7 7 (0^)+O(Ar 4 ), (6) 
where the fields r\ and 7 take the values: 



AtWK: 

6 * 4 

f=— 2,-1,1,2 



7 (±1) = 1 + v/6/3, 7 (±2) = 1 - V6/3 



(±1) =±W2(3- VI), ij(±2) = ±W2(3 + V6 



(7) 



Since ifj is invariant under a rotation in spin space, the above Hubbard-Stratonovitch decomposition conserves SU (2)- 
spin symmetry. Thus, the choice of the HS transformation for the Hubbard term will determine whether the algorithm 
is 5*C/(2)-spin invariant or not. 

We have carried out our simulations with the Projector QMC algorithm [0||. Within this approach, the ground 
state expectation value of an observable O is obtained with: 

(*o|0|*o> = j. (* T \e- eH Oe-° H \* T ) 

r\ /iTr l„ — 2ttffliTr \ V ' 



(*o|*o) (* T |e~ 2eH |* T ) 

The ground state l^o) is filtered out of a trial wave function, \^t), provided that (^oI^t) 7^ 0. We choose the trial 
wave function to be a spin singlet solution of the non interacting Hamiltonian (U=W=0). An explicit construction 
of such trial wave functions may be found in reference [g] . 

After Trotter decomposition of the imaginary time propagation and HS transformation of the two-body terms, one 
obtains: 

/; =^Pr(6,.x)(0)(e,x)+0(Ar 2 ), (9) 



(* T |e- 2eff |*T) 

where x denotes a configuration of HS fields, and (0)(9, x), corresponds to the value of the observable O the the HS 
fields x. At half-band filling particle-hole symmetry leads to positive values of Pr(0, x) which may thus be interpreted 
as a probability distribution and sampled with Monte Carlo methods. This statement is valid for both choices of the 
HS transformation of the Hubbard term (Q, |^). 

We now compare the SU(2) spin invariant algorithm based on Eq. (^|) to the SU(2) spin non-invariant algorithm 
based on Eq. (|l|). The SU(2) spin-invariant algorithm forces us to work with complex numbers. On the other hand, 
for many applications real numbers may be used for the SU (2) non-invariant code. For the comparison discussed 
below, and to keep the CPU time approximately constant, we have carried out twice as many sweeps for the SU(2) 
non-invariant code than for the SU (2) invariant code. We consider various observables at half-band filling. 

a) Magnetization. In the case of the SU (2)-invariant algorithm, and the above mentioned choices of the trial 
wave function one has: 

(d A cj A )(Q,x) = (d lCll )(Q,x). (10) 

Thus, the total magnetization, m z (q) — ^^e 1 ^ (n^ — n^i) is identical to zero for all values of the HS fields: 

(m,($))(e,aO=0. (11) 

On the other hand, the SU(2) non-invariant algorithm, equation ( |l0|) is not valid, and one obtains zero magnetization, 
only after summation over the HS fields. At q — (tt, it) = Q, L — 6, (n) — 1, U/t — 4 and W/t — one obtains after 
2 x 10 5 sweeps, (m z (Q)) — —0.16 ± 0.31. Thus for this trivial case, the advantage of the SU(2) invariant algorithm 
over the SU (2) non-invariant algorithm is infinite. 
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L = 6, W/t = 0.0, U/t = 4.0, (n) = 1, T = 
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FIG. 1. Imaginary time displaced spin-spin correlations, at Q — The 5*17(2) invariant code is based on the HS 

transformation of Eq. (S) and SU(2) non- invariant on Eq. (|l|). To keep the CPU time approximately constant between the 
two simulations, we have carried out twice as many sweeps for the SU(2) non-invariant code as for the SU(2) invariant code. 
Here, we have used periodic boundary conditions, $ = in Eqn. (H). 



b) Imaginary time displaced spin-spin correlations. Here we consider the quantity: S a (q, r)S a (— q, 0) where 
S a (q,T) — J2j e lq: 'e TH S^e~ rH , S" being the a-component of the spin operator on site j. The numerically stable 
computation of imaginary time displaced correlation functions within the Projector QMC algorithm is described in 
Ref. 0. In the SU(2)-invariant algorithm one obtains with the above mentioned trial wave function: 

(S a (q,T)S a (~q,0))(O,x) = (STtf,T)F(-f,0))(e,x). (12) 

Here, a, 7 run over the three components of the spin. In the case of the SU(2) non-invariant the above equation 
is valid only after summation over the HS fields, x. Fig. 1 plots the spin-spin correlations for (n) = 1, U/t = 4 
and W/t = on a 6 x 6 lattice. The half-filled Hubbard model is expected to show long-range antiferromagnetic 
order in the thermodynamic limit. Thus, jj(S(Q,t)S(—Q,0)) is should saturate to a finite value in the large L and 
t limits. Here, L is the linear size of the square lattice and N the number of sites. On a finite size lattice, a spin 
gap is expected, thus leading to an exponential decay in r of the considered quantity. If one compares the quantity 
(S(Q,t)S(—Q,0)) for both codes, (Fig. la and Fig. lc) it is clear that the SU(2) invariant code does much better. 
The large fluctuations in the case of the SU (2) non-invariant code may be traced back to the x and y-components 
of the spin-spin correlation function. In fact, considering only the z-component of the correlation function (Fig. lb) 
yields good results. The HS transformation of equation (|l|) conserves the z-component of the total spin, but not the 
other components. We now consider the z-component of the spin-spin correlations and compare both algorithms. We 
have plotted in Fig. Id the relative errors for both codes. 
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L = 8, W/t = 0.35, U/t = 2.0, (n) = 1, T = 0, $ = $ / 2 
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SU(2) Invariant code: 
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FIG. 2. Same as Fig. hi for the parameter set: W/t = 0.35, U/t = 2 and at half-band filling. Here we use antiperiodic 
boundary conditions set $ = $o/2. It is clear that the 5C/(2)-invariant code does substantially better at large values of r. 

Overall, the SU(2) spin-invariant code produces larger errors within the considered r-range. However, for the 
S't/(2)-invariant code the relative error, is to a first approximation independent of r. In contrast, the relative error 
for the SU(2) non-invariant code grows as a function of r. At large values of r we expect the the SU(2) invariant 
code to be more efficient. To confirm this statement, we consider the parameter set: (n) = 1, U/t = 2, W/t — 0.35, 
$ = $0/2 on a 8 x 8 lattice. The data is shown in Fig. 2. It is clear that the S'C/(2)-invariant code does much better 
at large values of rt. In Fig. 2, the z-component of the spin-spin correlations ( Fig. 2b) shows larger fluctuations 
than (S(Q,t)S(—Q,0)} ( Fig. 2c). Exactly the opposite is seen in Figs, lb and lc. 
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SU(2) invariant code 


SU(2) non-invariant code 


SU(2) non- invariant code 




<n z (r)IP>t(0)> 


(U x (t)U. xA (0)) 


(W(t)W-HO)} 





0.15055 ±0.00080 


0.15097 ±0.00073 


0.13837 ±0.01506 


0.2 


0.05902 ± 0.00058 


0.05920 ± 0.00051 


0.05845 ± 0.00387 


0.4 


0.02915 ± 0.00055 


0.02941 ± 0.00040 


0.02846 ± 0.00215 


0.6 


0.01474 ± 0.00046 


0.01547 ± 0.00032 


0.01373 ± 0.00254 


0.8 


0.00728 ± 0.00046 


0.00859 ± 0.00039 


0.00841 ± 0.00105 


1 


0.00397 ± 0.00038 


0.00480 ± 0.00053 


0.00394 ± 0.00104 


1.2 


0.00240 ± 0.00033 


0.00251 ± 0.00033 


0.00086 ± 0.00080 



TABLE I. Il-mode correlation functions at half-filling for U/t = 4, W/t = 0.35 on a 6 x 6 lattice. Data in the last two columns 
were obtained with the SU(2) non-invariant code. Within this algorithm, the x and y components of the Il-mode correlation 
function are identical. For the SU(2) spin-invariant algorithm, the results are independent on the considered components of 
the correlation function. To keep the CPU time approximately constant between the two simulations, we have carried out twice 
as many sweeps for the SU (2) non- invariant code as for the SU (2) invariant code. 
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c) E-modes. The II modes introduced in the 50(5) theory of the unification of antiferromagnetism and supercon- 
ductivity are defined by: II Q = J2p s s > d(p) c p + Q s { <ja<jV )s s ' c p>' ®- Here, a a corresponds to the Pauli spin matrix 
and for the case of <i-wave superconductivity, we consider g(p) = cos(p x ) — cos(p y ). The 5?7(2)-spm invariant code 
satisfies 

(H a (r)IL^(0))(Q,x) = (E7(r)IF'+(0)>(e |a! ). (13) 

In the case of the SU(2) non-invariant code, the above equation is valid for all values of 7 and a only after summation 
over the HS fields x. In the table, imaginary time II correlations are considered for both algorithms. For the SU(2) 
non-invariant code, substantial fluctuations in (n 2 (r)n z '^(0)) are observed at large values of r. On the other hand, 
the error-bars in (H x (t)I1 x '' (0)) are smaller. We again attribute the large fluctuations in the z component of the II 
mode correlations to fact that IF does not commute with the z-component of the total spin. In contrast, IP does 
commute with the z-component of the total spin. The SU (2)-invariant code shows good convergence, and the results 
agree with those obtained for the x-component of the II-mode within the SU{2) non-invariant algorithm. 

d) Single particle Green functions. We have not found any significant improvements in the fluctuations of 
the imaginary time single particle Green functions between the two algorithms. To be more precise, the fluctuation 
of X)<j ?( c ?<j( r ) c i ) are to a first approximation independent of the choice of HS transformation. Thus for both 

considered codes, the error-bars scale as C '/ \j 'N sweeps where C is a constant independent on the choice of the algorithm, 
and N sweeps denotes the number of sweeps carried out. 

In conclusion, we have considered an alternative HS transformation for Hubbard type models which conserves 
SU(2) spin symmetry. This algorithm requires the use of complex numbers and does not introduce a sign problem at 
half-band filling. In the case of SU(2) non-invariant algorithms based on Eq. (|l|), very large fluctuations can occur 
in the calculation of imaginary time displaced quantities of the form (A(t)A') when A does not commute with the 
z-component of the total spin. In the SU (2) invariant formulation, this pathology does not occur. We have compared 
the two algorithms in the worst case scenario where real numbers can be used for the SU (2) non-invariant algorithm. 
By keeping the CPU time constant for both codes we have shown that the long imaginary time behavior of the 
spin-spin correlations are obtained more efficiently with the SU(2) invariant code than with the SU(2) non-invariant 
code. Thus, the SU(2) invariant code is more efficient for the measure of spin gaps which may be extracted from the 
long imaginary time decay of the spin-spin correlations. More generally, it is an efficient code for the study of spin 
dynamics. The fluctuations of other quantities such as single particle green function, were invariant under the choice 
of the HS transformation. 
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